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This paper attempts to develop a philosophy for empirical time series analysis, involving 
the routine use of four data handling procedures (covariance estimation, .spectral estimation, 
autoregressivc model fitting and spectral estimation, and trend elimination and estimation) 
embodied in a computer program. 

The cross-spectral analysis of a pair of time scries, each consisting of 4000 observations, 
requires approximately 10 minutes on a 7090, including computation of covariances. Several 
examples of empirical time series analysis are given. 



1. Introduction 

The probabilistic theory of time scries is now 
extensively developed, built on the pioneering work 
of Wiener, Khintchine, Kolmogorov, Cramer, Loeve, 
and Karhunen in the 1930's and 1940's. The 
statistical theory also enjoys an extensive literature 
and a fair share of beautiful results. While the 
probabilistic theory can be pursued for the sake of 
its great beauty, it would be a mistake if the statisti- 
cal theory were to be developed only for its elegance. 
The ultimate aim of the statistical theory must be 
to provide data handling procedures for achieving 
the aim of time series analysis: synthesis of stochastic 
models which can be used to describe, and perhaps 
to control, the mechanisms generating each time 
series and relating various time series. 

For this reason, one may define a field, which 
may be called empirical time series analysis, with 
aims such as the following: 

(i) to develop efficient computer programs for the 
statistical treatment of empirical time series, paying 
especial attention to flexibility of input and output, 

(ii) to develop a philosophy, based on statistical 
theory, for judging and interpreting the statistical 
data reduction provided by the computer output, 

(hi) to provide experience in the small sample 
applicability of statistical procedures derived from 
asymptotic theory, 

(iv) to focus attention on the theoretical questions 
requiring further investigation, such as the problem 
of how to transform observed data to put them into 
a form where they satisfy the assumptions required 
to apply various statistical data reduction routines. 
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The goal of empirical time analysis is precisely that 
enunciated by Richard Hamming [1962] for numerical 
analysis: the seeking of insight rather than numbers. 
One should not expect empirical time series analysis 
to be a mere matter of grinding out answers. Rattier 
one desires to obtain a wealth of answers which 
can be tempered by imagination and judgement to 
achieve models for time series. 

My aim in this paper is to sketch an approach to 
empirical time series analysis whose basic attitude 
is: one should analyze the data in a number of ways, 
each corresponding essentially to a different possible 
model for the observed time series. Each analysis 
provides estimates of the parameters (or incompletely 
specified probability law characteristics) involved 
in the model. Comparing the analyses provides 
rough tests of hypotheses concerning which model 
provides a better fit to the data. 

2. Standard Models 

The standard model adopted for the analysis of 
a time series {X t {t), t=l, 2, . . .} is to write it as 
the sum 

X i (t)=rn i (f) + Y i (t) (2.1) 

of its mean value function 

mS)=E[XM] (2.2) 

and the residuals or fluctuations 

Y i (t)=X i (t)-m i (t). (2.3) 

In order to have the possibility of statistical 
inference from a single finite sample of a time series, 
one assumes that the residuals are covariance sta- 
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tionary [Parzen, 1962] in the sense that (for i, j=l, 2) 
there exist functions Rij(v) of integers v=0, ±1, ... 
such that 

Gov [X t (t), X j (t+v)]=E[Y i (t)Y j {t+v)]=R ij (v). 

(2.4) 

We call Rij(v) the co variance functions of the time 
series, and 

Pu(*) = ; P yLi/2 (2.5) 



{iJ„(0)#„(0) 



the correlation functions. 

The covariance functions Rij(v) always possess 
representations as Fourier-Stieltjes integrals. It is 
next assumed that they are Fourier transforms of 
spectral density functions /^(co): 



R«(v)=j' r e i9 *M<*)d<*, 



i= 



(2.6) 



The (auto) spectral density function /»(co) is a 
nonnegative even function of co: 



/„(-«)=/«(«)><), !C0|<7 



(2.7) 



The cross-spectral density function/^ (co) is in general 
complex valued; its real and imaginary parts are 
called respectively the co-spectral density, denoted 
Ctjfa), and quadrature spectral density, denoted g#(«). 
These functions possess the following properties 
(writing z to denote the complex conjugate of z) 



fij(-u) =/ (o>) =/,<(«) , 
2tf(— o))=— JtoM- 



(2i 



One can define for two jointly covariance station- 
ary time series X t {- ) a variety of spectral quantities, 
as follows: 



w tj (<»)-- 



ju (u)fjj (co) fii(u)fjj (co) 



(2.9) 



called the coherency between X t (-) and Xj(-) at 
frequency co; 

'^ u)J W (2 - 10) 

called the gain at frequency co of the predictor of 
X t (.) given X t (.); 



9j\i(u) z 






(2.11) 



called the gain at frequency co of the predictor of 
X } {.) given X t (-); 



^(co)=arc tan 






(2.12) 



called the phase difference between the two series 
at frequency co; 



fj\i(o>)=Mu){i-w ij (u)} : 



(2.13) 



called the error spectrum of the predictor of X 3 {-) 
given X t (.). 

In the sequel we shall discuss the question of 
estimating the mean value functions and spectral 
functions associated with time series. In terms of 
the spectrum, one can devise various physical 
mechanisms (especially filters) which might have 
generated the time series and which might be used 
to simulate them. 



3. Sample Covariance and Correlation 
Functions 

A basic step in empirical time series analysis is 
to form estimates of the covariance, correlation, and 
spectral functions. I believe it correct to say that 
these estimation problems do not as yet have 
generally accepted solutions. One of the aims of 
this paper is to stress the points at issue. 

The first point I desire to raise is that for the sake 
of developing a modular computer program we should 
adopt the following definitions for the sample 
covariance functions i?n(-) and R 2 2(-), and sample 
cross-co variance functions B 12 (-) and /? 2 i(-) : f° r 
i,j = l,2, 



i\ t=i 



»=0,1, 



For negative values of v we defne 

Rij(— v) = R ij (v). 



,N-1. 
(3.1) 

(3.2) 



For the sake of distinguishing the estimates from 
the population quantities they are estimating one 
should write R*{v) instead of Rij{v) to indicate 
that it is an estimate. For ease of writing we omit 
such asterisks. However, they should be inserted 
by the reader when discussing the properties of the 
estimates. 

Computing the sample covariance and cross- 
covariance functions will be one of the most time 
consuming aspects of a time series analysis of real 
data. It turns out that one can choose a convenient 
number V m&x , less than N, such that one need 
compute only 

Rijfy), 0=0,1, .• .,Fmax. 

One rarely chooses V m&x to be greater than 40 percent 
of the sample size N, and often it will be 0.25 N. 
Some considerations on how to choose Fm ax are 
discussed below. We call V max the covariance 
truncation point. / 
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From the sample covariance and cross-covariance 
functions, one forms sample correlation and cross- 
correlation functions, denned by 



* ij (v)=R i MI{R ii (0)R jj (0) 



U/2 



(3.3) 



for v = 0, 1, . . ., F max . As an example of the 
meaning of (3.3) note that 



Pu{v)-- 



J=\ 

N > 

t=l 



v=0, 1, . . ., N-l. 

(3.4) 



It should be noted that some writers would esti- 
mate Pij(v) by the ordinary correlation coefficient 
between the two series 

{X f (0,t=l,2,. ..,N-v}&nd 

{X,(t+v),t=l,2,...,N-v} (3.5) 

which is given by 

S (Xt(t)-Xt t J (XM+vj-Xf,.) 

t = \ 



{N- 



1 t=\ J 



defining 



Xf.o- 



1 



T,x&), 



X! 



N—v fei 

1 N 



1/2 

(3.6) 

(3.7) 
(3.8) 



The next two sections discuss the motivations for 
and implications of definitions (3.1) and (3.3), and 
why we oppose definition (3.6). 

4. Mean Subtraction 

It seems odd that we do not use in thejsample 
covariance function the deviations X(t)—X of the 
observations from the sample mean 

X4|>Y(0 

rather than X(t) itself. To explain our position on 
this point we introduce the notion of detrending. 

In analyzing a time series {X(f), t=l,2, . . . , N, } 
one usually adopts a model for X(-) of the form, 

X(t)=m(t)+e(t), 

where ra(-) is the mean value function, sometimes 
called the trend, and c(-) is the fluctuation or residual 



series about the trend assumed to be covariance 
stationary. 

One is interested in the trend for two purposes: 

(i) to estimate it, as an important part of the 
model one is fitting to the time series, 

(ii) to eliminate it in such a way as to obtain esti- 
mated values of the residual series e(-) so that its 
covariance and spectral structure can be estimated. 

The optimum procedure for estimating the trend 
depends on the co variance structure of the residuals. 
Consequently, in a sense, the problem of trend elim- 
ination needs to be solved before one can solve the 
problem of trend estimation. 

Detrending refers to subtracting an estimated mean 
value function m(t) from a time series X(t) to pro- 
duce a new time series X d (t)=X(t)—rn(t). We call 
X d (>) a detrended version of A'(-) about its mean 
value function. Detrending methods differ only in 
the way in which the estimated mean value function 
is formed. Two important ways in which one at- 
tempts to eliminate trend (or detrend) are mean 
detrending and linear detrending. 

Mean detrending. The estimated mean value func- 
tion is the arithmetic average of the time series 
observations: 



Mt)=if fj X(s)=X for t = l, 2, 



., N (4.1) 



The detrended time scries is 
X d (t)=X(t)-X; U 



=1,2, 



N 



(4.2) 



Linear detrending. The estimated mean value 
function is a regression line fitted to the time series 
by the method of least squares: 



N 



m(t)=X+b(t-t) for*=l, 2, 
where 



S* X (t)-tJ2 X(t) ib tX(t)-NiX 



J2t 2 -N(ty 

1 = 1 



The detrended time series is 
X a (t)=X(t)-[X+b(t-t)] 



(2V 2 -1)2V/12 



(4-3) 

(4.4) 
(4.5) 

(4.6) 



=1,2,... ,2V. 

(4.7) 



We now see why in forming sample covariances 
one should not necessarily subtract out the mean of 
a time series. We believe that one should assume 
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that the time series either have zero means or rep- 
resent residuals after detrending. Further, there 
will be occasions when automatically subtracting out 
the mean limits our ability to compute interesting 
estimates. We will often desire to consider a con- 
stant time series 



X(t)=l, all*, 



(4. 



which we do not desire automatically to replace by 
the zero time series. 

5. The Divisor Question: TV— v or N 

In the definitions given by (3.1) for R tj (v) we have 
chosen to divide by N (the number of terms in the 
series) rather than N—v (the number of terms being 
summed to form R tj (v)). Many researchers seem to 
prefer to divide by N—v on the grounds that it leads 
to an unbiased estimate of the true covariance 
Rij(v) (in the case of time series X t {-) with zero 
mean) : 



^[ 3 ^g^(O^(*+«)]=j^gB«(«')=-B«0 



:»). 



The suggestion that it is preferable to use the 
divisor N, rather than N—v, which is adopted in this 
work, is motivated by the following two considera- 
tions. Let us call R t j(v) with a divisor of N—v the 
unbiased estimate and Rtj(v) with a divisor of TV 
the biased estimate. 

One may show that 

(i) The unbiased estimate regarded as a function 
of v is not a positive definite function while the 
biased estimate is. This property is desirable for 
two reasons; first, because we are estimating a posi- 
tive definite function, and second, because it ulti- 
mately leads to nonnegative estimates of the spec- 
tral density function. 

(ii) One may show that in many cases the biased 
estimate has a smaller mean square error than the 
unbiased estimate. While this seems to be true in 
general, a rigorous proof has not as yet been found 
[see Schaerf , 1963, for proofs in various special cases]. 

In her Stanford Ph. D. thesis, Mirella Casini 
Schaerf [1963] shows that while the unbiased sample 
covariance function usually has a greater mean square 
error than the biased sample covariance function, 
neither seems to provide a really satisfactory esti- 
mate of the true covariance function. The sample 
covariance functions never damp out to zero, which 
we assume is the case for the true covariance func- 
tion. Schaerf examines various ways of modifying 
the sample covariance function to improve its prop- 
erties, with somewhat pessimistic results. 

Nevertheless it seems that when properly trans- 
formed the biased sample covariance function (in 
spite of failing to give a reasonable picture of the 
true covariance function) does yield estimates of the 
spectrum which provide a reasonable picture of it. 
It possesses the essential property which the sample 
covariance function must have (in order for its 



Fourier transform to be interpretable as a spectral 
density function), namely positive dejiniteness. 

The definition of the sample correlation function 
given by (3.6) is not positive definite. There seems 
to be no reason to use such an estimate. In what 
way is it relevant to our aims in time series analysis? 

6. Estimates of the Spectrum and Cross 
Spectrum 

The theory of estimation of the spectrum and 
cross spectrum is too extensive to be conveniently 
summarized here. We can only state the estimates 
which present theory seem to indicate should be 
formed as the first step to understanding the spec- 
trum (I have tried to summarize this theory in 
Parzen [1961 and 1964]). 

Three methods of computing spectra seem to be 
available for consideration: 

(i) the indirect or transjorm method, which Fourier 
transforms weighted co variances to estimate spectra, 

(ii) the direct or filter bank method, which esti- 
mates spectra as the variance and covariance of 
various filtered time series, 

(hi) the method of autoregressive spectral estima- 
tion. 

We discuss method (i) in this section and method 
(hi) in section 8. We do not discuss method (ii); 
see Ormsby [1961], Welch [1961], and Brillinger 
[1963]. 

Let Rij(v) denote either the sample cross covariance 
or the sample cross-correlation function. While in 
our opinion one should use the latter for ease of 
interpretation, the graphs one obtains have exactly 
the same shape in either case since they differ by 
constant factors (assuming (3.3) is used). 

As an estimate of the true cross-spectral density 
function / i; (co) one forms the estimate 



b 'k(j^B t M (6.1) 



Air \v\<M 



which depend on a choice of two quantities: 

(i) an integer M, called the truncation 'point of the 
spectral estimate (we usually choose several trunca- 
tion points in practice), 

(ii) a kernel k( . ) known as the lag window of the 
spectral estimate; its Fourier transform 

K(co)=^-\ e- {v »k(u)du (6.2) 

is the spectral window generator of the estimate while 

K M (o>)=±-T l e i >«k(^) (6.3) 

is called the spectral window of the estimate; it may 
be shown that approximately 



K M (u)=MK{Mw). 



(6.4) 
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There is a third choice to be made in forming the 
estimate /#(«), and this is the number of points on 
the interval to t at which it will be computed. 
We adopt the attitude thatjf^(co) should be computed 
for 



A T O V 



(6.5) 



where Q is an integer to be chosen. We call Q 
the spectral computation number. 

Spectral window generators. In this work, we 
use only the following spectral window generator: 



k(u) 



= l — Qu 2 +6u 3 , 

=0, u>\ 

=Jc(—u), u<0. 



0<?/<0.5 
0.5<?/<l 



(6.6) 



A theory of spectral window generators is developed 
in Parzen [1964]. It is shown that the kernel (G.6) 
always leads to nonnegative spectral estimates whose 
variance (when properly normalized for comparison) 
is slightly less than the variance of various other 
estimates considered. In particular a comparison 
is made with certain other kernels which have been 
widely used; see Blackmail and Tukey [1958], p. 98. 
One such kernel is 



k(u)=- (1+cos ttu), \u\ <1 



=0. otherwise. 



(6.7) 



Statistical sign ificance oj spectral estimates. Assu n i - 
ing normality of the observed time series, it may be 
shown that the variance of the estimated spectral 
density 



#(»)= 



I&^Km)^ 



(6.8) 



depends on the sample size N, the kernel k(-), the 
truncation point M, and the true spectral density 
fu(ca) as follows: 



Var[./;*(a,)]=^ 



«(»)4 f" **(«)** \ if 0<a>O 



if co=0 or 



(6.9) 



It should be noted that the equality in (6.9) is only 
approximate; from a rigorous point of view it 
should be written as a relation that holds in the 
limit. 

Note that the variance in (6.9) increases as the 
truncation point M increases. Therefore one is 



tempted to make M small. However the smaller 
M the larger is the bias of f* } (oi) as an estimate of 
/«(w). Thus arises the crux of the spectral estima- 
tion problem how to choose the truncation point 
so as to optimally compromise between a number of 
conflicting objectives. 

If instead of jfi(oi) one plots its logarithm log e 
f?i(u>) one achieves two advantages: 

(i) one magnifies the graph of the estimated 
density in regions where it is small, enabling one to 
more easily study the behavior of the spectrum at 
these frequencies, 

(ii) one achieves an estimate whose sampling 
properties are the same for all frequencies since by 
large sample statistical theory (6.9) implies 

Var [log,/;* («)]=^/ f~P(u)du\ if 0<co<tt, 



-jK-O*} 



if 00=0 Or C0 = 7T. 

(6.10) 



Whereas a confidence hand about the estimated 
spectral density has a width which varies with the 
height of the density, a confidence band about the 
logarithm of the density has a constant width: 



|lo<> 



iS(«)-logJ„(o))|<2^J"^(tt)rfu\ 



1/2 



(6.11) 

can be considered to be a confidence band for each 
a; in ()<co<7r of approximate confidence 95 percent. 



For the kernel (6 



•6),£ 



k 2 (u)du=0M. Therefore 



the right-hand side of (6.11), hereafter denoted A, 
is 0.33 for M/N=0M and 0.46 for M/N=0A0. 
The value of A for other ratios M/N is then easily 
approximated (thus, A is 0.46 VJ= 0.92 for M/N=* 
0.40). 

From (6.11) it follows that for each co in 0<a><7r, 
with approximate confidence 95 percent, 



,> .#(") -/«(*>) ^ .A 



1. 



(6.12) 



It is worth noting how quickly the limits in (6.12) 
increase as M/N increases; thus 



M/N 


e~ A -l 


e A -l 


0. 05 


-0. 3 


0.4 


. 10 


-. 4 


. 6 


. 20 


-. 5 


. 9 


.40 


-. 6 


1. 5 
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Only for M/N less than 10 percent is the percentage 
error of the estimated spectrum of reasonable size. 

The statistical significance of cross spectral esti- 
mates is discussed by Jenkins [1963], and Goodman 
[1963]. 

Choice of truncation point. It seems to be in- 
creasingly accepted among workers in statistical 
spectral analysis that the spectrum should be 
computed for several choices of truncation point M. 
As yet we do not have routine quantitative (or even 
qualitative) procedures for interpreting the spectra 
obtained from several choices of M. Such procedures 
are still under development. My experience has 
led me to feel that a good picture of the spectrum 
can be obtained by taking three truncation points 
Af 1 <Af 2 <M 3 . I choose M x to be an even number 
between 5 percent and 10 percent of the sample 
size N. I then choose M 2 = 2M 1 and M 8 =2M 2 . In 
section 9 we give examples of spectra computed with 
several truncation points. 

Choice oj spectral computation number. The spec- 
tral computation point Q has in the past frequently 
been chosen to be equal to the truncation point M. 
One can prove a sampling theorem to the effect that 
the estimated spectrum (which is a function of co, 
measured in cycles per unit of observation time, 
in the interval 0<co<0.5) can be recovered from its 
value at M equally spaced points. However, this 
recovery cannot necessarily be done by linear 
interpolation. If the graph of the estimated spec- 
trum is to be obtained by merely drawing line 
segments connecting the computed values, one needs 
to compute the spectrum at Q equispaced frequencies, 
where Q should be at least 2M and perhaps should 
be4M. 

If one uses 3 truncation points M l <CM 2 < CM h it 
has seemed reasonable to me to compute each 
spectrum at Q=M^ points. However, one should 
choose Q (approximately equal to M?) such that the 
frequencies which are multiples of tt/Q are of physical 
interest. For economic time series of monthly data 
we usually choose Q to be a multiple of 12. 

7. Computation Formulas for Cross-Spectral 
Estimates by the Transform Method 

The time series X x (-) and X 2 (-) whose cross spectra 
one is estimating are usually not the directly observed 
time series but rather the result of various detrending 
and filtering operations. Assuming that we are 
dealing with two time series ready for cross-spectral 
analysis, the following computations are performed. 

Step 1. Let N be the number of observations in 
the two series. Choose an integer T r max <A^ and 
compute, for v=0, 1, 2, ..., F max , the cross 
covariances 



Rn(v),R 22 (v),R 12 (v),R 21 (v) 



(7.1) 



and the cross correlations which we denote by the 
notation _____ 

Rii(v),R*(v),R ia (v),Rn(p). (7.2) 



Step 2. Choose three truncation points M h M 2} 
M 3 . For each truncation point, carry out steps 3 
to 5. 

Step 3. Choose an integer Q. Then for each 
frequency 



"=°'i ,2 i * 



(7.3) 



compute the quantities described in steps 4 and 5. 

Step 4. We prefer to compute normalized cross 
spectra (the transforms of cross correlations rather 
than of cross covariances) since they seem to be 
easier to graph and to compare. 

The normalized spectral density functions, co- 
spectrum, and quadrature spectrum 

/a («),/_(«), Cx_(«), Stafa) (7.4) 

at frequency to could be computed by the formulas 



fu(u>)-- 



= l\l El2 



(0)+gcos(tto)i v 



&>««}' 



(7.5) 



*»(«)=£{§ *»(°)+|§ COS M*(j£) 

[R 12 (v)+R 2l (v)] (7.6) 

2u(«)=^g sin M^(^)[4W-S 21 W] (7.7) 

where M denotes the truncation point. 

These spectral quantities are most efficiently com- 
puted not by using the explicit formulas above but 
by using an efficient procedure for evaluation of 
finite Fourier transforms due to Goertzel [I960]. 

Step 5. We next compute the following spectral 
quantities: 

Amplitude, 4i2^) = |/i 2 (co)H{k 12 (co)| 2 +|( Zl2 (co)| 2 } 1 / 2 

Phase, ^ 1 2(a)) = tan _1 {gi 2 (co)/ci2(co) } 

Coherence, T^i2(co) = |/i 2 (co)| 2 ^/ii(co)jf 2 2(co) 

Gain 1/2, #1/2 (co) =^4i 2 (co) -s-/aa(«) 

Gain 2/1, G 2/l (co)=A 12 (cc)-^M^)' 

Computation time. It may be worth explicitly 
noting two ways which we have found for reducing 
the computation time required. 

In computing covariances we call upon an external 
procedure for computing inner products which exe- 
cutes faster than the equivalent coding in program- 
ming language. Figure 7 a lists our co variance 
computation routine. 

In computing spectra and cross spectra we use 
the Goertzel procedure mentioned after (7.7). Fig- 
ure 7b lists this computation routine. 
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PROCEDURE COVARlANCE(NtMtLl»L2tX( ) $ R1()*R2( 
COMMENT THIS PROCEDURE COMPUTES THE AUTO AND 
TIONS* Rl(I) tR2( I) tCKI) AND CT(I)t FOR 1 = 
AT LAG ZERO IS STORED AT I=l» THE FUNCTION 
I=M+1. THE TIME SERIES ARE OF EQUAL LENGT 
THE ARRAY Y()» ONE BEGINNING AT LI* THE OT 
FUNCTIONS ARE NORMALIZED TO HAVE A VALUE 1 
CROSS CORRELATIONS ARE ALSO CONSISTENTLY N 
IZING FACTORS ARE D1#D2 AND 03. THE FUNCT 
ARRAYS Rl( )#R2( )tCI() AND CT ( ) TO ALLOW PO 
lNPROD(K f L»N»A( ) »B( )) IS AN EXTERNAL FUNCT 
'SUM*CU0t FOR I*(0»l f N-l)t SUM*SUM+A ( K+I ) • 
BEGIN 

INTEGER I..* 

Dl * INPROD( 

D2 = INPROD( 

D3 « SORTID1 

FOR KK * (It 



) *CI( ) tCTC )$01 tD2eD33$ 

CROSS CORRELATION FUNC- 
lt2t«««tM+l« THE FUNCTION* 

AT LAG M IS STORED AT 
H N AND BOTH ARE STORED IN 
HER AT L2« THE AUTO-CORR 

AT THE ORIGIN AND THE 
ORMALIZED. THE NORMAL- 
IONS ARE ADDED INTO THE 
OLING OF COVARIANCES» SUM» 
ION EQUIVALENT TO 
8(L+IH S 



,J.*.tK.«« *L... tM< 
LltLlfNtXC ) tX( ) )$ 
L2,L2,N,X( ) tX( ) )$ 
♦ D2) $ 
ItM+l) S 



,N« 



BEGIN 



END $ 



Rl (KK)=Kl(KK 
R2 (KK)=R2(KK 
CI (KK)=CI (KK 
CT (KK)=Cj (KK 

RETURN END $ 



)+INPROD(Ll »L1+KK-1#N-KK+1,X( ) tXf I I/D1S 
)+INPROD(L2»L2+KK-l f N-KK+ltX( ) f X( ) )/D2$ 
)+INPROD(ll,L2+KK-l.N-KK+l f X< ) »X( ) )/D3$ 
)+INPROD(L2 f Ll+KK-l,N-KK+l» X( ) tX( ) )/D3$ 



Figure 7a. Procedure covariance. 



The computation time T of the covariances of a 
pair of time series is approximately the product of N 
(the length of the time series) and M (the truncation 
point). Some typical values (on a 7090) are: 



N 


M 


T 


200 


80 


4 sec. 


1000 


200 


1 min. 


4000 


500 


8 min. 



The computation time T of a cross-spectral analy- 
sis is approximately the product of M (the truncation 
point) and Q (the spectrum computation number). 
Some typical values (on a 7090) are: 



M 


Q 


T 


40 




40 


1 sec. 


80 




80 


4 sec. 


80 




160 


8 sec. 


200 




200 


24 sec. 


500 




200 


1 min. 



8. Autoregressive Spectral Estimation 

Given auto regression coefficients a u . . ., a m one 

can form a new time series, denoted X(-), from a given 
time series X(-), by the formula 

X(t)=a l X(t-l) + . . . +a m X(t-m). (8.1) 

A 

We say that X(-) is obtained by autoregressive filter- 
ing from X(-). The residuals 



*{t)=X{t)-X(f) 



(8.2) 



are said to be autoregressive residuals. They are 

examined to determine how good a predictor X(t) 
is oiX(t). 

The autoregressive coefficients a i} . . ., a m may 
either be specified a priori or may be estimated from 
the data by a procedure known as stagewise auto- 
regressive estimation (see fig. 8a for the computation 
routine). In this case, one may be able to estimate 
the unknown spectral density function of the time 
series X(-) by a method called autoregressive spectral 
estimation. To describe this method we discuss some 
relations between transfer functions. 
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PROCEDURE TRANSFORM (M.N.Wl ) t Rl < ) tREO tROt ) tR2< )SF1 ( ) f F2 C ) »CO( ) »QUC ) ) * 
BEGIN 

COMMENT THIS PROCEDURE COMPUTES N+l POINTS OF TWO EST 
DENSITY FUNCTIONS AND OF THE CO-SPECTRUM AND QUADR 
FROM Rl() AND R2M WHICH ARE THE AUTO-CORRELATION 
AND RO< ) WHICH ARE THE EVEN ANO OOD PARTS OF THE C 
FUNCTION OF 1WO TIME SERIES. THE TRUNCATION POINT 
ING KERNAL USED IS W<)* THE SINES AND COSINES NEE 
RECURSIVELY, REFERENCES 1. )HAMMl NGtR.tf. t «NUMERICA 
SCIENTISTS AND ENG I NEERS • tMCGRAW-H I LL 1962.PAGES 7 
INTEGER Ift J ft <•••• I .... M...» N... * 
PI « 3»1415927 S PIV » 0.31830989 S 
C « N S CI • C3 « 01 ■ C0S<PI/C) $ 

D4 « C2 « SIN(PI/OS 06 * 2. CI J 
P1*0.5.R1(1) S P2*0.b.REm * 



IMATED SPECTRAL 
ATURE-SPECTRUM 
FUNCTIONS AND RE( ) 
ROSS-CORRELATION 

IS M. THE WEIGHT 
DEO ARE COMPUTED 
L METHODS FOR 
1-74 $ 



FOR I-(2tltM+l) % 
BEGIN 

A « W( I > S 
P1*PHR1( I ).AS 
END % 
Fl(l) « PI .PIV S 
C0( 1) ■ P2.PIV S 
FOR I « ( 1 .l.N)S 
BEGIN 

Ul) « U12 « U13 * 
FOR J ■ (M+lt-1.2 
BEGIN 
A - 
U31 



P?»P2+RE< D.A* 



P4=0.5.R2( 1) * 



P4«P4^R2 I I ) .A$ 



F2U ) 

0U( 1 ) 



U14 



P4.pfv 

o.o $ 



U21 



U22 



U23 



U24 « 0.0 * 



U3? 
U31 
U34 
Ull 
U13 
ENO 
Fl ( 1 + 1 > 
C0( 1 + 1 > 
0U< 1+1 > 
F2CI+1) 



W( J) $ 

« D6.U21-U11+R1 < J).A 

* D6.U22-U12*RE< J). A 

* D6.U23-U13+RO< J) .A 

* D6.U24-U14+R2C J).A 

* U21 S U21 » U31 

* U23 I U23 • U33 



U12 
Ul* 



-(D1.U21 ,- Ull 
«(D\.U2? - U12 
« D4.U23.PIV S 
MD1.U24 - U14 



RETURN 



01 ■ C1.C3-C2.D4 % 

C3 « 01 $ 
END % 
END $ 



♦ Rl ( 1) .0.5). PIV 

♦ RE( 1) .0.5) .PIV 

♦ R2< 1) .0.5). PIV 
D4 « D4.C1+C3.C2 S 
D6 ■ 2.D1 S 



U22 
U24 

% 
$ 



U22 
U24 



U32 

U34 



Figure 7b. Procedure transform. 



If the given time series X(-) actually was a white 
noise process, with spectral density function 



/x(«)=; 



<2ir, 



(8.3) 



then the time series e(-), defined by (8.2) would have 
as its spectral density 

jf.(«)=/x(«)-|A(«)|» (8.4) 

where A ((a) is the frequency transfer function 

4(a):=l— Otf-*"— . . . — a m <r rofw , (8.5) 



representing (up to a factor e iut ) the output of the 
operation defining e (•) when the input is e ifat . There- 



>peration 
fore 



fe(u) = — {(1— a x cos co- 
+ («i sin co+ . 



— a m cos mw) 2 



+a w sin mw) 2 



(8.6) 



Given autoregressive coefficients a X) . . ., a m we 
denote the right-hand side of (8.6) by TEW (co), and 
call it the w/wte noise transfer function of the scheme. 

We desire to normalize the white noise transfer 
function to have unit area (over the interval |co| <ir). 
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2 PROCEDURE SELECT* ROt K, L $COF(),IN(), NC )$ BEGIN 

2 COMMENT THIS PROCEDURE COMPUTES THE BEST FITTING AUTO-REGRESSIVE 

1 COEFFICIENTS FOR A TIME SERIES WITH SAMPLE CORRELATION FUNCTION R() 

2 OF TRUNCATION POINT K. L IS THE LENGTH OF THE TIME SERIES* USED FOR 
2 COMPUTING DEGREFS OF FREEDOM. THE NC SIGNIFICANT COEFFICIENTS 

2 ARE RETURNED IN COF ( ) WITH THEIR INDICES IN IN(). REFERENCES 1.) 

2 EFROYMSON ,M. A. ♦• MULTIPLE REGRESSION ANALYSIS' iIN 'MATHEMATICAL 

2 METHODS FOR DIGITAL COMPUTERS », ED I TED BY RALSTON, A. AND WILF, H.S., 

2 WILEY, 1962 S 

2 GLOBAL A( , ) $ 

2 INTEGER I.... J..., K..., L..., u. # ., N... S 

2 APRAY CB<72> $ 

2 ARRAY CHI2(5?) =(6.63, 9. 21, 11.3, 13.3, 15.1, 16.8, 18.5, ?0.1, 21.7, 

2 23.2, 24.7, 26.2, 27.7, 29.1, 30.6, 32.0, 33.4, 34.8, 36.2, 37.6, 

2 38.9, 40 # 3, 41.6, 43.0, 44.3, 45.6, 47.0, 48.3, 49.6, 50.9, 52.2, 

2 53.5, 54. 8» 56.1, 57.3, 58.6, 59.9, 61.2, 62.4, 63.7, 65.0, 66.2, 

2 67.5, 68.7, 70.0* 71.2, 72.4, 73.7, 74. 9, 76.2, 77.4, 78.6 ) S 

2 TOL = 0.001 $ N = K+l $ F2 = 6.635 5 KP = $ 

2 COMMENT PLACING THF COVARIANCE FUNCTION IN THE MATRIX $ 

2 FOR I = ( 1, 1 ,N) S 

2 BEGIN 

2 A ( I , I ) = 1 . $ 

1 FOP J*(lt1tl-1)$ A(I.J) = A(J,I) = R(I-J+1)$ 

2 END $ 
2 FIVE.. 

2 Fl = CHI2(K - KP ) $ PHI ■ L - KP -1$ RTOT = 0.0 $ 

2 VMIN = 2*30 $ VMAX = CO $ NMIN = NMAX = S 

2 FOR I=(1,1,K) S IF (A(I,I) GTR TOL )S 

2 BEGIN 

2 V = ( A( I ,N) ,A( N, I ) )/A( I , I ) $ 

2. IF V GTR 0.0 $ RTOT = RT0T + V $ COMMENT I NOT IN MODEL $ 

1 EITHFP IF (V LSS 0.0) S COMMENT I IS IN MODEL $ 

2 BEGIN 

2 CB( I )=A( I »N) $ 

2 IF (ABS(V) LSS ABS(VMIN) )$ (VMIN= VS NMIN = IS )$ 

2 END S 

2 OR IF (V GTR VMAX)$ (VMAX= VS NMAX=I$ CB(I)= 0.0S)S 

2 OTHERWISE S CB( I ) = 0.0 S 

2 ENDS 

2 EITHER IF ABS(VMIN) LSS F2 . A ( N ,N ) /PH I S (K1 = NMIN $ KP=KP~1S )$ 

2 OR IF RTOT GTR Fl . ( A ( N ,N ) /PHI ) $ (Kl = NMAX $ KP=KP+1 S )% 

2 OTHERWISES GO TO ALM S 

2 COMMFNT WE COMPUTE THF NEW MATRIX. IF Kl = NMIN THEN WE ARE DELETING A 

2 VARIABLE. IF Kl = NMAX WE ARF ADDING A VARIABLE. IN EITHER CASF THE 

2 COMPUTATIONS ARE THE SAME S 

2 PIVOT = 1.0/ A(K1,K1) S 

2 FOR J=(1,1,N)$ A(KltJ) = A(K1 , J) .PIVOT S COMMENT CHANGE PIVOTROWS 

2 FOR I=( 1,1 »N)$ IF I NEO Kl $ 

2 BEGIN 

2 TEMP = A( I ,K1 ) S 

2 FOR J«(l,ltN)$ A(I,J)= A(I,J)- TEMP.AIK1 ,J )$ 

2 A( I ,K1 )=-TEMP. PIVOT $ 

2 ENDS 

2 A(K1,K1 > ■ PIVOT S 

2 GO TO FIVE $ 

2 ALM.. 

2 NC = S 

2. FOR I=(K,-1,1)$ 

2 IF (CB(I) NEQ 0.0)$ (NC = NC + 1S COF ( NC ) =CB ( I ) $ I N ( NC ) =K- 1 + 1 $ ) $ 

2 IF NC GTR 20 $ NC = 20 S 

Z RETURN END $ 



Figure 8a. Procedure select. 
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One may verify that 

f * TRW (o)dct=l+al+ . . . +a%. (8.7) 

Therefore the normalized white noise transfer function, 
denoted by TRWN («), is given by 



TRWN («) = 



TRW (») 



l+a?+...+a* 



(8.8) 



Note that this transfer function is normalized to 
have unit area over the interval — t<oo<t. 

If X(-) is in fact an autoregressive scheme satis- 
fying the model 

X(t)=a 1 X(t-l)+ . . . +a m X(t-m)+e(t), (8.9) 

where e(-) is white noise, then from (8.4) its spectral 
density function is given by 

/x(«)=y.(«)H-|ii(«)l» 

= [27r{ (1— a x cos co— . . . —a m cos mco) 2 

+ (&! sin co+ . . . -\-a m sin mco) 2 }] -1 . (8.10) 

We denote the right-hand side of (8.10) by 
TRAR (a>), and call it the autoregressive transfer junc- 
tion of the scheme of coefficients a h . . ., a m . 
Note that 



TRAR (co)-{4x TRW (co) J- 1 . 



(8.11) 



We have not been able as yet to find a formula 
(convenient for computation) for the integral 

TRAR (co)rfco. Since we prefer to compute nor- 



/: 



malized transfer functions, in practice we approxi- 
mately evaluate this integral by a crude method of 
numerical integration and thus compute the normal- 
ized autoregressive transfer function, denoted by 
TRARN (co) and given by 



TRARN ( w ) = 



4ttTRW (co); 



f'{^rTRW (co)}" 1 ^ 



(8.12) 



We can describe the method of autoregressive 
spectral estimation. It consists of (i) determining 
autoregressive coefficients a lf . . . a m by stage wise 
autoregressive estimation, (ii) estimating the spectral 
density function / e (co) of the autoregressive residuals, 
(hi) if y e (co) is approximately the spectrum of white 
noise, we take the normalized autoregressive transfer 
function, denned by (8.12), as the estimated spec- 
trum of the original time series X(t). 

The statistical theory of stagewise autoregressive 
model fitting is considered by Schaerf [1963]. 



9. Some Examples of Empirical Time Series 
Analysis 

One of the major aims of empirical time series 
analysis is to provide a modern solution to what 
classically has been called "the search for hidden 
periodicities. " Classical approaches to this problem 
have been discredited because they seemed to pro- 
vide evidence for the existence of "spurious cycles." 
One way to avoid seeing cycles in data where they 
are not present is to compare estimated spectra from 
this data with the estimates formed under the same 
sample size and truncation point from data of known 
properties. 

As an example of the above method, an analysis 
of an economic time series will be carried through. 
The radio case would be essentially the same, apart 
perhaps from changes in length and time scale in the 
series, and the economic material has already been 
analyzed by our method. 

Given a series of length 180 (as many economic time series 
are) one might choose a minimum truncation point of 16; the 
other truncation points would then be 32 and 64. 

To obtain some idea of the resolution of our spectral win- 
dows, it is instructive to first perform a spectral analysis of a 
constant time series 



X l (t) = l, t=l, 2, 



N. 



(9.1) 



No mean subtraction is performed. The sample covariance 
function is given by 

The resulting estimated spectral density functions are given 
in figure 9a. The curves corresponding to different truncation 
points are easily distinguished since at zero frequency the 
curves increase in amplitude with increasing truncation point. 
The spectral density function corresponding to M=16 has a 
local maximum of 0.175 which might be construed as a 
"spurious" cycle. 

Series of stock price indexes have been extensively analyzed 
(see Granger and Morgenstern [1963]). We consider a 
monthly series of such an index, denoted X 2 (t), for the 180 
months in 1948-1962, taking logarithms and mean detrending. 
The estimated spectra, given in figure 9b, are exactly as in 
figure 9a, except for a little additional power at high fre- 
quencies. 

Similarly we estimated the spectra of the linear series 



X 3 (t) = t, 



= 1,2, 



180 



(9.3) 



with no detrending. The estimated spectra were again 
exactly as in figure 9a. 

The conclusion to be drawn is that for both the stock price 
index series and the linear series, the spectrum consists essen- 
tially of a single line (or spike) at zero frequency; a more pre- 
cise statement of this assertion might perhaps be made using 
techniques of mixed spectral analysis such as those currently 
being investigated by George Hext [1964]. 

To characterize more clearly the mechanism generating the 
stock price index series and the linear series, other than to 
say they have a spectrum concentrated at zero frequency, one 
fits these series by autoregressive schemes using stagewise 
autoregressive estimation. For the stock price index series 
one finds the scheme 



X 2 (O = 0.986X 2 (*-l) + e 2 (O 



(9.4) 
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Figure 9a. Log of estimated spectral density functions, constant series. 
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Figure 9b. Log of estimated spectral density functions, stock price index 48-62, mean detrend. 
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Figure 9c. Lo(/ o/ estimated spectral density functions, stock price stageauto. 




Figure 9d. Log of estimated spectral density functions, Fed. cash from public 48-62 mean detrend. 

948 





30.0- 


















! 1 1 - 1 ! ••• ' • i '■'■ ■'■■ i FMr;TH ikn ■ 














; : '-:j ■ j . . i i 




.:.!.:. 










J' 1. ! i- : ::| 'j 1 .'_ I i ''■- ! H9 












i i .!'.•' 


• -i i - '---1 — --f--— 




mrJiB 


: .::;; 


1 






! | ■ : ; ' 1 




.:. ']■ : 


-! - 






1 i ■• ; ■ i i 


\ ■ 1 ■ j : j .' 


■ " j i- ; 1 ! ' - i : " 






::!- 








■..: ■• j ; J 


i J : i--.-S--i 


j ; 














! : • ' . 1 


. ', . | .':' : 






-r :: - 










■ ' : I i ! i : \ ' : ■ ' 




: j_. - ; - - 






\ j : ■ • 




l. o - 




1 i ! • i ; ! : ' : ' : 


j ! 1 : 














~:- 








I 


: \ 1 - : ! 






















i . ; 1 : ; 1 


! • ■■ 1 


■ : ] 












: . : 








: ! ' : ■ ! . 


' ! . 




"i 1 ' • ' 


! - ' 


















~-~- :i w^ ; i 


■;.:.■ 






1 -""'•' 






-— 


— 






:1- 














. ■ | 


j 






; | 






. :.:.L..- ■_ 


VT^ 










?m 'r 4- - 




;\.l ;. 


. 1 








§§ 








: ... \ 


fc^H^^ 


. ^iV- 




•S ,1- : 1 ; ! 


' 


1 }>* 


»N : '-'il 




n 


}.2B- 




mi 


-■; j : : - ■ 


\J::fr 


MR 


fjfjgi 














iK 




4JD ■ 


: i : ! : 






I m 






$mM !v/i 

.: | .: I :. i . : ! ■ 




















: : 








■ 




•!■:!■! i 










































■ ! ! : 




; 




. 1 : ■• 




■ 








1 ! ■ i i ■ 










i i • 








3. 01 - 






! " ! 


1 : 
























n : 


















j 






i 


"'-^ i • 




















• ! • 




: 1 : : 








! 1 i ■ 






















: ' : 












i ;-. ] 










■ | ' ! ! i 1 


, ! [r ^ -f- 












i 1 : 1 -■ ■ 1 : ' 




. ! . . 


















i • i 1 - 






























001 - 


_i^_ 










;■ m 


—j 


m 


1 : i : i i i 


i 




s 





FREQUENCY IN CYCLES PER UNIT TIME 

Figure 9c. Log of estimated spectral density functions, Fed. cash stageauto. 



while for the linear series 



^ 3 (O = 0.992A r 3 (i-l) + € 3 (O. 



(9.5) 



We then estimate the spectra of the residual time scries 
e 2 (t) and e 3 (0- 

The estimated spectral density functions of the stock price 
index residuals arc given in figure 9c; they are well within the 
confidence limits of white noise spectra. The true spectrum 
of white noise on the logarithmic plot used in the figure is a 
horizontal line at 



Loi 



«.(iooo£). 



:5.07 



(9.6) 



From the foregoing analysis we have essentially obtained 
the often found random walk model for the logarithms of 
stock market prices X 2 (t): 



X 2 (t)=X 2 (t-l) + 62(t) 



(9.7) 



where 62(0 is a white noise series. It would be interesting to 
investigate whether the models (9.4) and (9.7) are significantly 
different. 

The estimated spectral density functions of the linear series 
residuals are again exactly as in figure 9a. Examination of 
the printed residuals e 3 (Q shows that they are given by a 
linear series 



X 3 (*)=l+(<-D0.00831, t=l, 



(9.8) 



Performing a stagewise autoregression on this series disclosed 
that it approximately satisfies the model € 3 (£)=0.994e 3 (2— 1); 
one would conclude that the original series X s (t) satisfies 
X z {t) — 2X z (t— 1)4- X :i (t-2)=0 and thus is a straight line. 



Our last example is intended to illustrate how autore 
gressive spectral estimation is used to check the results of 
the transform method of spectral estimation. Let X±(t) 
denote a monthly series of logarithms of cash payments re- 
ceived by the Federal Government from the public. Its 
spectra are graphed in figure 9d. There are prominent peaks 
at w = and co = 0.33. There are questionable peaks at 
0.083, 0.167, 0.25, 0.42 (for the interpretation of these fre- 
quencies in terms of periods see the table below). 

Fitting an autoregressive scheme to X±(t) by stagewise 
autoregressive estimation one finds the scheme 



X 4 (0 = O.732X 4 (*- 



■ 3) + O.825X4 (t- 12) 

-0.611X 4 (Z- 



■15)4-€ 4 (0. (9.9) 



The spectra of the residuals e 4 (0, given in figure 9e, are not 
quite white; there is a slight predominance of power at low 
frequency which should be investigated. Nevertheless, it 
is white enough for us to feel justified in taking the normalized 
autoregressive transfer function as an estimate of the spectrum 
of X 4 (t). This transfer function is given in the 1/TRANSFER 
column of figure 9f. All the peaks previously found are 
still present, with their relative size clearly indicated. The 
model given by (9.9) lends itself readily to prediction. The 
relative mean square prediction error is defined to be 

Var[e 4 (Q] 
Var [X 4 (t)] 

These variances are routinely computed in our program. Their 
ratio turns out to be about 0.1; the predictor is thus rather 
reliable. 

Associating frequencies and periods. Consider a time series 
observed at monthly intervals. The following table shows 
the frequencies corresponding to several important periods. 
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AUTOREGRESSION COEFFICIENTS USED IN DETRENDING ARE 


.7322 X(T- 3) 








.8252 X(T-12) 








-.6108 X(T-15) 


REQ. 


TRANSFER 


1/TRANSFER 


LOC(TRANS) 


.0000 


.0002 


3.8218 


-8.6503 


.0104 


.0030 


.2214 


-5.3020 


.0208 


.0185 


• C362 


-3.9906 


.0312 


.0553 


• 0121 


-2.8949 


.0417 


. 1025 


.0065 


-2.2777 


.0521 


. 1271 


.C053 


-2.062C 


.0625 


. 1018 


.0066 


-2.2851 


.0729 


.0401 


.0167 


-3.2163 


.0833 


.0028 


.2403 


-5.8840 


.0937 


.0571 


.0117 


-2.8636 


.1042 


.2165 


.0031 


-1.5299 


.1146 


.4123 


.0016 


-.8860 


.1250 


.5285 


.0013 


-.6577 


.1354 


.4841 


.0014 


-.7254 


.1458 


.2998 


.C022 


-1.2046 


.1563 


.094 1 


.0071 


-2.3629 


.1667 


.0054 


. 124L 


-5.2228 


.1771 


.0941 


.0071 


-2.3629 


.1875 


.2998 


.0022 


-1.2046 


.1979 


.4841 


.0014 


-.7254 


.2083 


.5285 


.0013 


-.6377 


.2188 


.4123 


.0016 


-.8860 


.2292 


.2165 


.0031 


-1.5299 


.2396 


.0571 


.0117 


-2.8636 


.2500 


.0028 


.2403 


-5.8840 


.2604 


.0401 


.0167 


-3.2163 


.2708 


.1018 


.C066 


-2.2851 


.2813 


.1271 


.0053 


-2.0625 


.2917 


. 1025 


.0065 


-2.2777 


.3021 


.0553 


.0121 


-2.8949 


.3125 


.0185 


.0362 


-3.9906 


.3229 


.0030 


.2214 


-5.8020 


.3333 


.0002 


3.8218 


-8.6503 


.3438 


.0030 


.2214 


-5.8020 


.3542 


.0185 


.0362 


-3.9906 


.3646 


.0553 


.0121 


-2.8949 


.3750 


.1025 


.0065 


-2.2777 


.3854 


.1271 


.0053 


-2.0625 


.3958 


.1018 


.0066 


-2.2851 


.4063 


.0401 


.0167 


-3.2163 


.4167 


.0028 


.2403 


-5.8839 


.4271 


.0571 


.0117 


-2.8636 


.4375 


.2165 


.0031 


-1.5299 


.4479 


.4123 


.0016 


-.8860 


.4583 


.5285 


.0013 


-.6377 


.4688 


.4841 


.0014 


-.7254 


.4792 


.2998 


.0022 


-1.2046 


.4896 


.0941 


.0071 


-2.3629 


.5000 


.0054 


.1241 


-5.2228 



LOG( 1/TRANS) 

1.3407 
-1.5076 
-3.3190 
-4.4147 
-5.0319 
-5.2471 
-5.0245 
-4.0933 
-1.4257 
-4.4461 
-5.7797 
-6.4236 
-6.6720 
-6.5843 
-6. 1050 
-4.9467 
-2.0869 
-4.9467 
-6. 1050 
-6.5843 
-6.6720 
-6.4236 
-5.7797 
-4.4461 
-1.4257 
-4.0933 
-5.0245 
-5.2471 
-5.0319 
-4.4147 
-3.3190 
-1.5076 

1. 3407 
-1.5077 
-3.3190 
-4.4147 
-5.0319 
-5.2471 
-5.0245 
-4.0933 
-1.4257 
-4.4461 
-5.7797 
-6.4236 
-6.6720 
-6.5843 
-6. 105O 
-4.9467 
-2.0869 



Figure 9f . Transfer functions. 



950 



Period 9, length 


Frequency = 1/0 


of a cycle 
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2 


. 500 


3 
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6 
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12 
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It is a pleasure to thank Howard Taylor who 
contributed the bulk of the computer program that 
we use for empirical time series analysis. I desire 
to express my appreciation for the great intelligence, 
ingenuity, persistence, and patience that he brought 
to bear on this work. 
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